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ABSTRACT 

A  time -dependent  numerical  scheme  is  used  to  solve  for  the 
chemical  nonequilibrium  profile  behind  a  normal  shock  wave  in  air. 
The  steady  state  equations  are  also  integrated  using  a  fourth  order 
Runge-Kutta  technique  and  a  comparison  is  made  with  the  time- 
dependent  calculation.  The  results  agree  within  one  percent  except 
in  the  region  close  to  the  shock.  In  this  region  the  profiles  differ 
because  the  time -dependent  technique  allows  calc-  lation  through  the 
shock  (which  is  several  mesh  points  in  width)  and  as  a  result  some 
dissociation  occurs. 
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I.  INTRODUCTION 


A  good  deal  of  effort  has  been  devoted  to  studying  gas  dynamic 

problems,  including  thj  effects  of  chemical  reactions.  Earlier  work 

generally  involved  systems  with  a  single  dissociation -recombination 

reaction.  A  survey  of  flow  studies  of  that  class  has  been  given  by  Li 

Mure  recent  investigations  have  involved  the  study  of  gas  dynamic  flows 

2 

with  a  series  of  coupled  chemical  reactions.  The  work  of  Marrone  is 
an  example  of  calculations  of  this  type.  These  results  have  all  been 
obtained  by  solving  the  steady  state  equations. 

In  this  paper  we  use  a  time  dependent  technique  to  solve  for  the 
nonequilibrium  region  behind  a  normal  shock  wave  moving  in  air.  We  do 
not  consider  the  discontinuities  to  be  internal  moving  boundaries  but  rather 
obtain  their  motion  from  the  solution  of  the  differential  equations.  Hence 
the  methods  used  to  integrate  the  differential  equations  yield  weak  solutions 
in  the  space  time  domain. 

As  initial  conditions  for  this  problem,  we  choose  a  discontinuous 
function  consisting  of  two  constant  equilibrium  states.  The  fluid  particles 
at  equilibrium  ahead  of  the  shock  move  through  the  shock  into  the  equili¬ 
brium  state  behind.  The  concentrations  of  the  various  species  within  the 
particle  are  no  longer  in  equilibrium  with  the  new  temperature  and  pressure 
and  chemical  reactions  begin  to  occur.  Since  the  velocity  of  the  fluid  behind 
the  shock  is  less  than  the  shock  velocity,  the  distance  between  the  two 
original  equilibrium  states  increases  until  sufficient  time  has  elapsed  to 
bring  the  entire  system  to  a  steady  state,  relative  to  the  shock.  The  gas 
is  assumed  to  be  thermally  perfect  and  in  local  thermodynamic  equilibrium. 
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In  a  previous  paper  ^  one  of  us  presented  results  for  a  calculation  of 
this  type  for  a  Lighthill  gas  using  a  first-order  scheme.  In  this  paper  we 
use  a  scheme  which  is  more  accurate  and  do  the  calculation  for  a  six  re¬ 
action  model  for  air. 


II.  CHEMICAL  MODEL 


The  i**1  of  r  reactions  involving  species  j  may  be  written  as 


s  K{{  S 

Zai  .M.  -  7  b.  M. 

J  v -  **  J 

j*l  %  j=l 


(1  -  1, . . . .  r/ 


where  the  a;.  andbi.  are  the  stoichiometric  coefficients  of  speciee  j.  M- 
J  J  J 

is  the  chemical  formula  for  species  j  and  Kf^  and  are  the  forward  and 
backward  reaction  rate  constants  for  the  i  reaction. 

The  molar  rate  of  production  of  species  j  due  to  reaction  i  is  given 
by 

°,  .*(b.  .-a. .)  Kf. 
ij  '  ij 


8  '  a  s  k. 

7 T  Caaia-^7T  C  Dla 

q=  i 


■1  n  — 


a=I 

where  Kc  is  the  equilibrium  constant  based  on  molar  concentrations  and 
the  Ca  ’  s  are  the  molar  concentrations. 

The  total  rate  of  production  of  species  j  is  the  sum  of  the  tates  of 
production  of  species  j  over  all  r  reactions. 

VI  % 

i=I 

For  each  atom  in  the  mixture  there  will  be  an  equation  of  the  form 


(1) 


(2) 


L  j  lJ 


.th 


where  ^  is  the  number  of  atoms  of  the  j  element  per  molecule  of  the 


(3) 
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i^1  species.  This  relationship  expresses  the  conservation  of  the  j^1 

element  and  allows  computation  of  the  0 .  for  elements  in  terms  of 

J 

the  for  molecules. 

The  following  six  reactions  are  the  assumed  chemical  model 


for  air: 


o2+Mj=£20+m 

N2+M^=i2N+M 


NO+O: 


n2+0- 


iM+O+M 


i°2+N 

1NO+N 


n2+o2-^bA2no  (f) 

The  species  M  in  the  first  three  reactions  may  be  02>  0,  N,  NO 
or  Ar.  In  general,  there  is  a  different,  forward  rate  constant  associated  with 
each  M.  The  forward  rate  and  equilibrium  constants  used  have  the  following 


form: 


K  =Af  T^f  exp  (  -  -m-  ) 


Kc=AcT*Cexp  (  -  -^) 


Values  for  the  constants  A^,  Ac>  q^,  qc»  E0^,  and  E0  are  listed  in  Table  I. 

The  first  three  reactions  involve  dissociation  by  two  body  collisions 
and  lead  to  the  formation  of  atoms  by  the  direct  dissociation  of  the  molecules. 
The  reverse  reactions  lead  to  molecular  formation  through  atom  recombination 
by  three  body  collisions.  Reactions  (d),  (e),  and  (f)  are  bimolecular  exchange 
reactions. 
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in.  DIFFERENTIAL  EQUATIONS 


The  equations  to  be  solved  may  be  written  in  the  following  vector 


form: 


w  =f  +  S 
t  x 

where  w  is  a  vector  function  of  x  and  t. 


(4) 


and 


0 


The  momentum  m  and  total  energy  E  are  defined  per  unit  volume. 

m2 

The  total  energy  is  the  sum  of  the  kinetic  energy  and  the  internal 

energy  pe. 

2 

m 

EsPe  + 

The  total  man  density  p  ia  Riven  by: 


■  %  wici 


where  W.  ia  the  molecular  weight  of  specie*  j. 
J 

For  our  equation  of  state  we  use 

P  =  CTR  T 


where 


'T=  £  °j 

j=l 
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A  specification  of  m,  E  and  the  Cj  allows  determination  of  p 
from  Eq.  (6),  T  from  Eq.  (11),  and  p  from  Eq.  (7).  Hence  we  have 
a  closed  system. 

IV.  INITIAL  CONDITIONS  AND  DIFFERENCE  SCHEMES 


The  initial  condition  consists  of  a  uniformly  moving  discontinuity 
separating  two  equilibrium  states.  The  upstream  eq*-  librium  state  and 
propagation  velocity  of  the  discontinuity  are  specified.  The  Rankine- 
Hugoniot  relationships  are  then  used  to  obtain  the  corresponding  downstream 
state: 

PlVP2U2  <12> 

Pl+piu12--p2+p2u22  (13) 

hl+J«  u^  =  h^^2  (14) 


In  the  above  equations  let  subscripts  1  and  2  indicate  the  upstream 
and  downstream  states  respectively,  u^  and  u^  are  the  velocities  relative 
to  the  discontinuity.  Severed  iterations  are  required  to  obtain  state  2. 

Using  Eqs.  (12)  and  (13)  we  may  write 

—  =  1  +  (1  -  —  )  (15) 

P  2  Y  jMj  Pj 


7 


where  Y  is  the  ratio  of  the  specific  heat  at  constant  pressure  to  the  specific 
heat  at  constant  volume  and 


Mr 


W 

Vi 


is  the  relative  upstream  Mach  number. 

From  the  Eq.  of  state  (7)  we  obtain 


TZ  _  P  I  W2 

Tl’Pl  °2  *1 


(16) 


where 


C.W. 

i  l 


For  the  first  step  of  the  iteration  procedure  choose  a  value 
of  PzfPl  anrt  calculate  P2/P  i  from  (15).  Solve  for  T2/T  j  using  equation 
(16)  with  the  assumption  W2=Wj.  With  the  values  of  T£  and  thus 
obtained,  calculate  the  corresponding  equilibrium  composition.  Since 
there  are  five  possible  products  being  considered  as  a  result  of  chemical 
reaction  between  0 ^  and  there  will  be  five  unknown  quantities;  namely, 
the  resulting  concentration  of  each  species.  The  five  equations  used  to 
obtain  these  unknowns  are  the  two  atom  conservation  equations  for  oxygen 
and  nitrogen  and  the  three  equations  relating  the  equilibrium  constants 
for  the  stoichiometric  dissociation  of  0^  and  ^  and  formation  of  NO  to  the 
concentrations  of  each  species. 

Using  the  calculated  concentrations  obtain  a  new  value  of 
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From  equation  (16)  obtain  a  new  value  of  T2.  Repitition  of  this  iteration 
will  yield  a  T^.  and  composition  satisfying  equations  (15)  and  (16). 

To  determine  whether  T2  satisfies  the  energy  equation  (14),  write 


h2  =  ^  + 


<AL,2] 

p  2 


hj  is  calculated  from  equation  (9)  and  hence  is  known.  We  may  also 

tl 


write  h2  as 


h_=  — ~ 

2  p2 

which,  after  substituting  equation  (10)  for  FL  yields 


h2=^ 


£l  <Ci>2AiT2  +  £  <Ci>2BiTZ  4=1  <Ci>2H?i 


This  equation  is  solved  for  If  T2  from  equation  (17)  is  not  equal  to  the 

one  found  previously  assume  a  new  value  of  ^2^1  31141  repeat  the  procedure. 

The  initial  value  for  18  calculated  assuming  frozen  flow  and 

fixed  Y  .  This  value  is  then  increased  until  the  above  conditions  are  met. 

Time -dependent  difference  methods  that  may  be  used  for  the  solution 

of  equations  (4)  have  been  discussed  by  one  of  us  in  Ref.  5  .  The  two-step 

method  used  in  this  calculation  is  given  below.  In  the  first  step  a  first 

order  scheme  is  used  to  compete  the  intermediate  values  at  twice  t  +  At. 

4 

We  denote  these  intermediate  values  by  a  bar.  A  second  order  accurate 

scheme  is  used  to  compute  the  final  values.  The  overall  scheme  has 

second  order  accuracy.  The  second  step  averages  the  f  differences  at 

t  and  t  +  &t  so  that  values  of  w,  f,  and  S  are  all  centered  at  the  point 
A  t 


(n,  t  + 


). 


-  t+At  .  t  t.,,.,t  -t  .  A  At  .  et  ie,t  . 
V  I,  %(wn+l+wn'iX<fn±  (Sn+Sn±l> 


- 1+  At 

7 

n 


w  =  Uw  .  .  +  w  , 

n+  1  n- 1 


)  + 


Crf;.i> +  t 


41  (S*  ,+S*  ,) 

'  n+ 1  n-  1 


The  bars  signify  itermediat".  values. 


(17) 
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The  final  values  are  computed  using 


t+At  _  t  \ 

W  =  W  +  — »r- 

n  n  Z 


+  ¥-  ^*  +  3 

Z  n  n 


_t+^t  t+At' 

L(f,  ,  -r  ,)  +  f  . ,  -7  . 

*'n+l  n-1'  n+  J5  n-lj 


t+A  t 


where 


,  -  At 

X"33T 


(18) 


The  initial  conditions  given  below  are  for  Mj  =  10.03  with  upstream 
conditions  representative  of  150,  000  feet  above  sea  level. 


p/l 

a 

9.67 

VP1 

= 

126. 75 

VT 1 

= 

11.  91 

*02 

= 

7.  2169X10-2 

xo 

a 

1.  9291X10” 1 

-1 

= 

6.8441X10 

*N 

= 

1.  7476X10”4 

*NO 

a 

4.  2199X10”2 

XAr 

= 

8.  131 1X10 

TH »  X.  are  the  mole  fractions,  i.  e. ,  number  of  moles  of  species 
i  per  mole  of  mixture. 

C. 

v  _  i 


Since  the  Initial  conditions  are  equilibrium  states  all  dependent 
variables  remained  constant  ahead  of  the  shock  and  in  the  region  behind  the 
relaxation  zone. 


ID 


5  6 

It  is  well-known  '  that  the  finite  difference  solution  to  Equations 

(4),  with  S  =  0,  for  numerical  schemes  of  the  type  employed  here  have  an 

overshoot  in  flow  variablos  immediately  behind  the  shock.  Inclusion  of 

chemical  reactions  does  not  eliminate  this  overshoot  since  their  effect  is 

negligible  through  the  shock.  The  overshoot  in  pressure,  in  particular, 

yields  a  value  which  exceeds  that  in  the  equilibrium  region  downstream. 

A  fluid  particle  experiencing  this  numerical  overshoot  will  relax  along  a 

dil.ereni  thermodynamic  path  and  hence  yield  a  different  nonequilibrium 
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profile.  We  'hose  to  eliminate  the  overshoot  by  using  a  first  order 
scheme  in  the  region  close  to  ihe  shock.  The  additional  numerical  vis¬ 
cosity  of  this  scheme  smooths  out  the  overshoot  but  lowers  the  accuracy 
of  the  calculation  in  this  region.  The  enfre  nonecuilibrium  region  extends 
over  two-hundred  mesh  points.  The  first  order  scheme  is  used  at  ten. 

v.  discussion  OF  RESULTS 

Results  presented  in  figures  I  through  VIII  correspond  to  a  Mach 
number  of  10.03  and  upstream  conditions  representative  of  *50,  000 
feet  above  sea  level.  The  calculation  was  performed  using  a  time  step 
of  1.  175X10  ^seconds  and  reuults  are  plotted  at  time  step  six-hundred. 
Comparison  with  time  step  five-hundred  indicated  no  change  in  the  profile. 
Hence,  these  results  are  for  a  steady  state  relative  to  the  shock.  To 
check  convergence  a  run  was  made  for  twice  the  number  of  time  steps  using 
half  the  value  of  At.  Results  were  the  same  to  within  one  percent.  A  further 
check  on  the  solution  was  afforded  using  a  fourth  order  Runge-Kutta  scheme 
to  integrate  the  steady  state  equations.  Jump  conditions  based  on  frozen 
omposition  were  used  at  the  initial  condition  for  this  calculation.  Compari¬ 
son  between  the  steady-state  and  time -dependent  techniques  showed  the 
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to  be  in  agreement  to  within  a  percent  except  in  a  region  close  to  the 
shock.  The  results  differ  in  this  region  because  the  time -dependent 
technique  allows  calculation  through  the  shock  (which  is  several  mesh 
points  in  width)  and  as  a  result  some  dissociation  occurs.  Thus,  the 
jump  conditions  obtained  are  not  based  on  frozen  composition.  The 
jump  conditions  based  on  frozen  composition  and  those  obtained  using 
the  time -dependent  method  are  listed  below  for  comparison. 


Frozen 

Calculated 

p^p, 

=  121.0892 

119.  7179 

VT1 

=  17.  3466 

16.3469 

VP1 

=  6.  9806 

7. 1774 

CM 

>f 

=  .21 

. 174726 

xo 

=  Q 

.039803 

*N2 

=  .781 

. 754165 

*N 

=  0 

.000126 

=  0 

.022359 

XAr 

=  .009 

.008821 

Figure  1  shows  the  effect  of  the  chemical  system  on  the  temperature. 
Immediately  behind  the  shock  the  temperature  is  450C°K.  It  then  decreases 
monotonically  through  the  relaxation  region  to  its  equilibrium  value  of  about 
3300°K.  The  total  density  increases  about  30%  from  its  value  immediately 
behind  the  shock  (Figure  2  )  whereas  the  pressure  remains  practically 
constant  after  a  slight  initial  rise  (Figure  3). 

Figure  4  shows  the  decrease  in  mole  fraction  of  as  tl»e  distance 
from  the  shock  increases  and  Figure  5  shows  the  corresponding  formation 
of  the  atomic  species  0.  The  molecular  species  N ,,  decreases  monotonically 
behind  the  shock  (Figure  6)  but  N  and  NO  have  maxima  (Figures  7  and  8). 
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The  NO  maximum  has  been  observed  experimentally  The 

reason  for  the  maximum  in  the  NO  and  N  mole  fraction  follow;)  from  an 
analysis  of  the  reaction  rates  through  the  non  equilibrium  region.  Because 
of  the  initial  abundance  of  0o  and  N9  reaction  (f)  is  at  first  a  major  produces 
of  NO.  Since  reaction  (a)  is  relatively  fast  Og  decomposes  into  0  atoms 
quickly.  As  0  appears  reaction  (e)  begins  production  of  NO  and  N.  Reaction 
(d)  initially  proceeds  to  the  left  making  this  reaction  a  major  contributor  to 
the  amount  of  NO.  Farther  downstream  of  the  shock  reactions  (e)  and  (f) 
reverse  and  proceed  to  the  left  until  equilibrium  is  obtained.  Reaction  (d) 
also  reverses  and  proceeds  to  the  right  as  equilibrium  is  approached. 

When  this  occurs  NO  and  N  begin  to  diminish.  Hence  we  get  maxima  in 
their  mole  fractions  at  approximately  the  same  point  behind  the  shock. 

For  inviscid  flows  with  chemistry  the  Runge-Kutta  integration  of 
the  steady  state  equations  is  faster  and  slightly  more  accurate.  Tha  advan¬ 
tage  of  the  time -dependent  method  lies  in  its  ease  of  extension  to  viscous, 
heat  conducting,  and  radiating  flows.  With  dissipative  effects  there  is  no 
numerical  overshoot  and  the  calculation  may  be  done  without  the  introduction 
of  the  lower  order  scheme  in  the  region  of  the  shock.  We  hope  to  report 
soon  on  a  calculation  of  this  type. 
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TA8I£  I 


Reaction 

FORWARD  RATE 

Catalytic  Body 

AND  EQUILIBRIUM  CONSTANTS 

Af  9f 

Ac 

9c 

E„c  (°K) 

(a) 

Ar,  NO,  N 

3.6 

(19)  -1.0 

5.95  (4) 

1.2  (3) 

-0.5 

5.95  (4) 

n2 

4.8 

(20)  -1.5 

5.95  (4) 

1.2  (3) 

-0.5 

5.95  (4) 

°2 

1.9 

(21)  -1.5 

5.95  (4) 

1.2  (3) 

-0.5 

5.95  (4) 

0 

6.4 

(23)  -2.0 

5.95  (4) 

1.2  (3) 

-0.5 

5.95  (4) 

(b) 

Ar,  0^,  NO,  0 

1,9 

(17)  -0.5 

1.13  (5) 

18.0 

0.0 

1.13  (5) 

n2 

4.8 

(17)  -0.5 

1.13  (5) 

18.0 

0.0 

1.13  (5) 

N 

4.1 

(22)  -1.5 

1.13  (5) 

18.0 

0.0 

1.13  (5) 

(c) 

Ar,  N2,  02 

3.9 

(20)  -1.5 

7.55  (4) 

4.0 

O.C 

7.55  (4) 

NO,  0,  N 

7.9 

(21)  -1.5 

7.55  (4) 

4.0 

0.0 

7.55  (4) 

(d) 

— 

3.2 

(9)  +1.0 

1.97  (4) 

3.3  (-3) 

+0.5 

1.60  (4) 

(e) 

— 

7.0 

(13)  0.0 

3.80  (4) 

4.5 

0.0 

3.75  (4) 

(f) 

— 

4.6 

(24)  -2.5 

6.46  (4) 

1.35  (3) 

-0.5 

2.15  (4) 

Note:  (N)  =  10N 
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TABLE  II 


COEFFICIENT  FOR  ENTHALPY  CURVE  FIT 


Species 

A^  (erg /mole /°K2  ) 

(erg/mo le/°K) 

H^i  (erg/mole) 

°2 

7.3947250136(+3) 

3.254044833l(+8) 

0 

0 

2 .4291027  537 (+2 ) 

2.1092339308(+8) 

2 .4677232 (+12  ) 

N2 

8. 23444622 42 (+3) 

3. 107 52 07 522 (+8) 

0 

N 

2 . 586549470l(+3 ) 

2 . 004883069l(+8 ) 

4. 7 078368 (-112  ) 

NO 

7.8127087774(+3) 

3.1875545777(+8) 

0.8985977(+12) 

Ar 

0 

2.0786723C00(+8) 

0 
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ALTITUDE  *150,000  FT. 


TEMPERATURE  VS  POSITION 


FIG.  2  TOTAL  DENSITY  VS.  POSITION 
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FIG.  3  PRESSURE  VS.  POSITION 


Vi  ^  *  10.03 
ALTITUDE  =  150,000  FT 
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FIG.  4  MOLE  FRACTION  0,  VS.  POSITION 


oo  =  10.03 
ITUDE>  150,000  FT 


FIG.  6  MOLE  FRACTION  N„  VS.  POSITION 


M  eo  s  10.03 

ALTITUDE  *150,000  FT. 


POSITION  (Cm) 

FIG  7  MOLE  FRACTION  N  VS.  POSITION 


M  co  *10.03 

ALTITUDE  *150,000  FT. 


Si 

o 


Si  2 

§  g 


FIG.  8  MOLE  FRACTION  NO 
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